Three-dimensional solitons and vortices in dipolar Bose-Einstein condensates 
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Three-dimensional solitary and vortex structures in Bose-Einstein condensates are studied in 
the framework of Gross-Pitaevskii model including the simultaneous action of local cubic-quintic 
nonlinearity and nonlocal dipole-dipole interactions. Nonlocal interactions are shown to change 
significantly the formation threshold and the numbers of atoms confined into the coherent structures. 
An appearance of robust high-order (m — 2) three-dimensional vortices is revealed. 
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INTRODUCTION 



Recent progress in experimental and theoretical stud- 
ies of Bose-Einstein condensates (BECs) has opened up 
a new opportunity to investigate nonlinear interactions 
of atomic matter waves [HQ. By applying an external 
magnetic field and using the Feshbach resonance, the s- 
wave scattering length can be tuned, thus it is possible 
to explore extreme regimes of interaction from strongly 
repulsive to strongly attractive. The resulting nonlinear 
evolution of matter waves gives the possibility to observe 
the nonlinear effects such as atomic self-focusing and for- 
mation of solitons. 

Dark [3] and bright solitons [3, 5, 6] have already been 
created in BECs. In the absence of an external trap, 
such 3D structures are always unstable: they either col- 
lapse, if number of atoms exceeds some critical value, or 
spread out in the opposite case. During an implosion 
of BEC, the atom density becomes high, and repulsive 
three-particle interactions come into play and should be 
taken into account, which gives rise to the additional local 
quintic nonlinear term in corresponding Gross-Pitaevskii 
equation (GPE) InRef. jlH 3D non-spinning soli- 

tons and vortex solitons have been studied in the frame- 
work of nonlinear cubic-quintic (CQ) Scrodinger equa- 
tion with the application to the bulk nonlinear optical 
media. The competition between the cubic attractive 
and quintic repulsive terms is able to arrest the collapse 
and to stabilize 3D solitons. In contrast to non-spinning 
solitons, vortex solitons can be unstable with respect to 
azimuthal perturbations which brake 3D vortex into sev- 
eral filaments. In conservative CQ media, single-charged 
3D vortices can be stabilized, while higher-order vortices 
remain unstable [ll| . Very recently, stable double-charge 
vortex solitons were found in the framework of the com- 
plex Ginzburg-Landau equation [l2[. As the matter of 
fact, dissipative solitary structures may be more com- 
pact and more robust, thus it may be easier to find 
stable localized vortices in dissipative systems than in 
their conservative counterparts. In the present paper we 
present the first example of stable 3D spinning higher- 
order (m > 1) solitons in a conservative nonlinear sys- 
tem. The stabilization is achieved due to the action of 
the nonlocal nonlinearity associated with long-range in- 



terparticle interactions in dipolar BEC. 

Nonlocal nonlinear media response naturally appears 
in a wide variety of physical systems such as plasmas 
HE E0, HI, Bose- Einstein condensates BEC 0, op- 
tical media (25[, liquid crystals [23|, Hij], an d s °ft mat- 
ter [27l ]. In the 2D case, nonlocal nonlinear interactions 
were shown to suppress an azimuthal instability and to 
stabilize vortex solitons [l5|, , [IK • It was found re- 

cently [28[ that spinning 3D solitons are unstable against 
splitting into a set of stable fundamental solitons in the 
medium with nonlocal thermal nonlinearity. 

In degenerate dipolar gases, the nonlocal nonlinear- 
ity arises from long-range, partially attractive, and 
anisotropic dipole-dipole (DD) interactions. The theoret- 
ical investigations have shown that the stability of dipo- 
lar gases crucially depends on the trap geometry @,l8j] . 
However, it is not always necessary to include an external 
trap since self-action of BEC can be enough to provide 
the conditions for formation of stable localized coherent 
structures. In the present paper we address the following 
question: how the nonlocal DD interactions affect the 
parameters and stability properties of the self-consisted 
solitons and vortex solitons formed in the trap- free BEC. 



II. MODEL EQUATIONS 

We consider BEC of atoms with electric dipole d (the 
model is equally valid for magnetic dipoles) oriented in 
the z direction by a sufficiently strong external field, and 
that hence interact via a dipole-dipole potential 



V d (r)=g d (l-3cos 2 tf)/r 3 



(1) 



where g& = ad 2 /A-keq, Sq is the vacuum permittivity, d 
is the angle formed by the vector joining the interacting 
particles and the dipole direction. The parameter a can 
be tuned (see e.g. [26[) in the range —1/2 < a < 1. 

In the mean field approximation, a dipolar BEC at 
sufficiently low temperatures is described by a GPE with 
nonlocal nonlinearity: 



+g d * [ V d {v-v')\^{v',t)\ 2 d z v' = Q 



(2) 
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where g = A-kTi 2 cl/M, a is the s- wave scattering length. 
In the following we consider attractive two-particle in- 
teraction (a < 0) and repulsive three-particle interaction 
(9k < 0). Note that we account here only for the con- 
servative part of three-particle interaction. Thus, GPE 
([2]) conserves the norm (the number of atoms in the con- 
densed state): 



TV 



Vr, 



(3) 



energy: 



E = 



A_iv*i> + >|<-I 9if |*i« 



- J7,il«f ©J <i»r, 9 = J Vi(r - r')|*(r')| 2 <iV, (4) 

momentum and angular momentum. 

In the next section, the general properties of stationary 
solitons and vortices are studied by analytical variational 
method and numerically. 



III. STATIONARY 3D SOLITONS AND 
VORTEX SOLITONS 

Stationary solutions of the Eq. ([2j) have the form 
^(r, t) = ^(r) exp(iAt) and obey the dimensionless equa- 
tion: 



Xijj + Ail; + i/; |^| - ip |^| + Cif) 6 = 0, 



V d (k)F 



(5) 



(6) 



where F denotes the Fourier transformation, Vd(k) — 
^-(3k 2 /k 2 — 1) is the Fourier transform of dipole-dipole 
interaction potential. Applying the following rescal- 
ing transformations: r dim iess = (2Mg 2 /h 2 g K ) 1/2 r p h y s, 
^dimiess = \qk I g | x ^ 2 V^phys ? the number of free parameters 
is reduced to two, namely C = gd/g, and A = —p\gK\/g 2 
(/i being the chemical potential). 

In order to study the general properties of station- 
ary 3D solitary and vortex solutions analytically, we em- 
ployed the variational analysis with the following ansatz: 



*(r,t) 



h I — 



|m| 



exp 



P 2 



(7) 



obtain the energy functional (|4|) at fixed number of atoms 
TV as the function of two variational parameters a p and 
a z . These calculations were performed analytically, but 
the resulting expressions are too cumbersome to be pre- 
sented here. Stationary solutions correspond to the sta- 
tionary points of the energy functional E, i.e. the pa- 
rameters a p and a z are found from the set of equations: 
dE/da z =0, dE/da p = at fixed TV. The dependen- 
cies TV (A) obtained by variational method below are com- 
pared with the numerical results. 

To find the stationary solutions of Eq. ([5]) numeri- 
cally, the stabilized iteration procedure similar to that 
proposed by Petviashvili (see, e.g. [29]) was implemented 
on a fully 3D grid with the resolutions up to 128 3 . Fig- 
ure [I] shows the isosurfaces of constant \ip\ 2 for vortex 
solitons m = 2 at different C. The change in sign of 
nonlocality parameter C leads to flattening (C > 0) or 
to elongation (C < 0) of the distribution of atoms along 
z-direction compared to the case (7 = 0. Below we shall 
concentrate mainly at the case C > 0. 

In the Fig. [2] the numbers of atoms versus the pa- 
rameter A for solitons (m = 0) and vortices (m = 1,2) 
are shown. The numerically found states are plotted by 
circles, while the dashed lines present analytical depen- 
dencies calculated by the variational method. It is seen 
that the variational predictions are in a good agreement 
with the results of numerical simulations. The diver- 
gence between analytical and numerical curves increases 
as the number of atoms TV confined into the structure 
goes up, and this can be explained by the observation, 
that with increase of A, soliton and vortex solutions de- 
velop more and more pronounced plateau and their pro- 
files then strongly deviate from the trial Gaussian profile 
used in the variational approach. In the Fig. O numbers 
of atoms are plotted versus the parameter A for different 
strengths of non-locality for the vortices m = 1 (depen- 
dencies TV (A) for solitons m = and vortices m = 2 show 
very much similar behavior). One can see that increased 
nonlocality leads to the significant decrease of the struc- 
ture's formation threshold and to the substantial elonga- 
tion of the range of accessible parameter A. 



(B) 





where p — ^/ x 2 + y 2 , tp is azimuthal angle. Inte- 
ger m is the topological charge, m = corresponds 
to non-spinning solitons, m > - to vortex soli- 
tons. Two variational parameters a p (t) and a z (t) de- 
scribe radii of the structure in the directions across 
and along the external field, and the amplitude h(t) = 
TV 1 / 2 7r _3 / 4 a~ 1 (t) [a 2 (t)m!] -1 ^ 2 is found from normaliza- 
tion condition (j3]). Using the trial function one can 



FIG. 1: Numerically found stationary vortex solutions (iso- 
surfaces of constant |V>| 2 ) for m = 2 and different values of the 
nonlocality parameter C: (A) C = 0, A = 0.1, (B) C = 0.3, 
A = 0.45, (C) C = -0.3, A = 0.08. 



3 



600- 



500- 



400- 



300- 



200- 



100 ■ 




equation which describes evolution of small perturbation: 
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FIG. 2: Number of particles N vs A for different stationary 
solutions at C — 0.3. Numerical results are shown in cir- 
cles, dashed lines are for variational predictions. The arrows 
indicate stability thresholds for vortices. 
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FIG. 3: Number of particles N vs parameter A for vortices 
with m — 1 with different C. Numerical results. 



IV. STABILITY AND DYNAMICS 



We start investigation of the stability of vortices with 
analysis of small perturbations applied to the stationary 
solution: 



where |e(r,t)| << |^|- Linearizing the GPE ([2]) in 
vicinity of stationary solution one gets the nonstationary 



(8) 



rif 

^-Ae + Ae+(2|Vfe + V>V) 
at 

-(2^V + 3|Vfe)H 2 + C(50 x/> + 6e) = 0, 
where is given by Eq. ([6]) and 

SG = F- 1 \v d (k)F {^e* + ^e} . 



We have solved Eq. ([8]) numerically by means of split- 
step technic. The azimuthal instability growth rate was 
calculated as follows: 



2At \ v(t) J ' 



where At is the time step, and v(t) = (e|e) is the norm 
of the perturbation. When the perturbation grows ex- 
ponentially, j(t) saturates at some value 7, which was 
taken as the maximum growth rate. 
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FIG. 4: 




and m = 2, C=0.3 

The soliton structures having zero topological charge 
are azimuthally stable, and their stability region coin- 
cides with that described by the Vakhitov-Kolokolov cri- 
terion (dN/dX > 0) [19]. The maximum growth rates 
obtained by means of above approach are given in Fig. [H 
for vortices m = 1 and m = 2 at C = 0.3. 

Note that our analysis is formulated in very general 
form, it gives the maximum growth rate at given A but 
not the growth rate of specific unstable azimuthal mode. 
This explains the knees in the dependencies 7(A) (see Fig. 
2j) which correspond to the intersections of growth rates 
7l(A) for different azimuthal numbers L of perturbations. 

As it is known (see, e.g., [30]), the competition between 
cubic and quintic nonlinear terms leads to the formation 
of solitary structures having a kind of plateau on the top 
when number of atoms N is high enough. Further in- 
crease of N leads to saturation of the amplitude and to 
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elongation of the plateau. In the previous works pjj[l4|, 
stabilization of the 2D vortices was associated with the 
appearance of such plateau. It is remarkable that in the 
BECs described by the GP model with additional non- 
local nonlinearity, stabilization of vortices occurs before 
the plateau is formed. At given nonlocality parameter C, 
the growth rates vanish at some threshold value of A due 
to the action of nonlocal dipole-dipole interactions. For 
instance, at C = 0.3, the single-charged vortices get the 
stability window at A > 0.26 and double-charged vortices 
(m = 2) - at A > 0.42, these threshold values are marked 
by the arrows in Fig. [2j As for vortices with m > 2, no 
stabilization was observed. 

The evolution and the dynamical stability of solitary 
and vortex structures was simulated numerically using 
the well-known split-step Fourier method [20[ , where the 
nonlocal DD integral term was calculated in the spectral 
space. We used perturbed solutions found by the station- 
ary solver as an initial condition for the nonstationary 
problem. For the structures with plateau, Petviashvili's 
approach may fail and is unable to reproduce the AT (A) 
diagram completely in the high-energy region. In this 
case, we used the appropriate initial condition with the 
parameters extracted from the variational analysis as an 
input for the non-stationary solver to check an existence 
and stability of solitons and vortices for values of N where 
our relaxative solver is inefficient. 
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FIG. 5: Decay of unstable vortex (m = 2, C = 0.3, A = 
0.25): isosurfaces of constant particle density \ip\ 2 are given 
for different times. 
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FIG. 6: Decay of unstable vortex (m = 2, C = —0.3, A = 
0.08): isosurfaces of constant particle density \ip\ 2 are given 
for different times. 

Direct numerical simulations confirmed the predictions 
obtained by the linear stability analysis. Solitons were 
found to be stable everywhere in the region dN/dX > 



0, stability thresholds for vortices coincide with those 
predicted by the linear stability analysis. In the unstable 
region, destruction of a vortex may occur in the different 
ways, depending on the number of atoms confined into 
the structure and on the applied perturbation. Typical 
unstable evolution snapshots are plotted in Fig. [5] and 
Fig. [6] for different signs of nonlocality parameter C. 

V. CONCLUSIONS 

In the present paper we study the influence of nonlocal 
dipole-dipole interparticle interaction on 3D solitons and 
vortices in BEC described by GPE with local contact at- 
tractive two-particle and repulsive three-particle interac- 
tion. Nonlocal anisotropic dipole-dipole interaction leads 
to quantitative, as well as to qualitative changes in scales, 
energies and stability properties of supported solitary 
structures. The formation threshold (minimum number 
of atoms, which is needed to form a structure) goes down 
compared to the case of absence of DD interactions, and 
possible range of chemical potentials significantly elon- 
gates. While soliton structures having zero topological 
charge are stable if dN/dX > 0, vortices may exhibit az- 
imuthal instability. We proved the stabilization of 3D 
vortices with m = 1 and m = 2 and revealed their sta- 
bilization thresholds. All higher order (m > 2) vortices 
were found to be unstable. Up to our knowledge, this is 
a first example of 3D stable vortex with the topological 
charge larger than one in a conservative system. 
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